Analizzare i dataset relativi all’inquinamento
Federico Reali
24/05/2017
1 Inquinamento a San Andreas
Iniziamo scaricando il dataset dal link dove ho isolato (e pulito) i dati raccolti dall’agenzia americana per la tutela dell’ambiente, relativi alle rilevazioni orarie della concentrazione delle PM10 nell’aria nella città di San Andreas, CA.
Questo dataset è una parte di un dataset molto più grande che si può trovare al link.
Una volta salvato il file nella nostra cartella data passiamo a leggere il suo contenuto.
Possiamo iniziare a prendere confidenza con il dataset leggendo alcune informazioni essenziali, come la dimensione, il nome delle variabili, oppure ne possiamo esplorare la struttura, visualizzare le prime o le ultime righe.
dim(PM10dataSA)[1] 2581 28Vediamo subito che il dataset contiene più di 2mila osservazioni per 26 variabili. Accediamo i nomi delle variabili per vedere se riusciamo a capirne il significato.
names(PM10dataSA) [1] "X.2" "X.1" "X"
[4] "State.Code" "County.Code" "Site.Num"
[7] "Parameter.Code" "POC" "Latitude"
[10] "Longitude" "Datum" "Parameter.Name"
[13] "Date.Local" "Time.Local" "Date.GMT"
[16] "Time.GMT" "Sample.Measurement" "Units.of.Measure"
[19] "MDL" "Uncertainty" "Qualifier"
[22] "Method.Type" "Method.Code" "Method.Name"
[25] "State.Name" "County.Name" "Date.of.Last.Change"
[28] "DateTime.Local" Proseguiamo visualizzando la struttura ed alcune righe dei dati.
str(PM10dataSA)'data.frame': 2581 obs. of 28 variables:
$ X.2 : int 1 2 3 4 5 6 7 8 9 10 ...
$ X.1 : int 1 2 3 4 5 6 7 8 9 10 ...
$ X : int 308088 308089 308090 308091 308092 308093 308094 308095 308096 308097 ...
$ State.Code : int 6 6 6 6 6 6 6 6 6 6 ...
$ County.Code : int 9 9 9 9 9 9 9 9 9 9 ...
$ Site.Num : int 1 1 1 1 1 1 1 1 1 1 ...
$ Parameter.Code : int 81102 81102 81102 81102 81102 81102 81102 81102 81102 81102 ...
$ POC : int 3 3 3 3 3 3 3 3 3 3 ...
$ Latitude : num 38.2 38.2 38.2 38.2 38.2 ...
$ Longitude : num -121 -121 -121 -121 -121 ...
$ Datum : Factor w/ 1 level "WGS84": 1 1 1 1 1 1 1 1 1 1 ...
$ Parameter.Name : Factor w/ 1 level "PM10 Total 0-10um STP": 1 1 1 1 1 1 1 1 1 1 ...
$ Date.Local : Factor w/ 120 levels "2016-01-01","2016-01-02",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Time.Local : Factor w/ 24 levels "00:00","01:00",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Date.GMT : Factor w/ 122 levels "2016-01-01","2016-01-02",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Time.GMT : Factor w/ 24 levels "00:00","01:00",..: 9 10 11 12 13 14 15 16 17 18 ...
$ Sample.Measurement : int 25 48 26 14 7 14 11 16 17 11 ...
$ Units.of.Measure : Factor w/ 1 level "Micrograms/cubic meter (25 C)": 1 1 1 1 1 1 1 1 1 1 ...
$ MDL : int 4 4 4 4 4 4 4 4 4 4 ...
$ Uncertainty : logi NA NA NA NA NA NA ...
$ Qualifier : logi NA NA NA NA NA NA ...
$ Method.Type : Factor w/ 1 level "FEM": 1 1 1 1 1 1 1 1 1 1 ...
$ Method.Code : int 122 122 122 122 122 122 122 122 122 122 ...
$ Method.Name : Factor w/ 1 level "INSTRUMENT MET ONE 4 MODELS - BETA ATTENUATION": 1 1 1 1 1 1 1 1 1 1 ...
$ State.Name : Factor w/ 1 level "California": 1 1 1 1 1 1 1 1 1 1 ...
$ County.Name : Factor w/ 1 level "Calaveras": 1 1 1 1 1 1 1 1 1 1 ...
$ Date.of.Last.Change: Factor w/ 4 levels "2016-05-06","2016-05-12",..: 1 1 1 1 1 1 1 1 1 1 ...
$ DateTime.Local : Factor w/ 2580 levels "2016-01-01 00:00:00",..: 1 2 3 4 5 6 7 8 9 10 ...head(PM10dataSA)tail(PM10dataSA)Vediamo che i dati contengono indicazioni come coordinate geografiche e lo stato dove sono stati raccolti i dati. Nel caso del dataset originale, queste informazioni sono essenziali per individuare il luogo di raccolta dati, mentre nel nostro caso potrebbero essere omesse, dopo aver verificato che siano consistenti in tutto il dataset.
Siamo interessati a capire l’andamento delle PM10 nel nostro dataset durante il tempo in cui sono stati raccolti i dati. Per prima cosa convertiamo i dati riguardanti le date nel giusto formato, cioè in Date. Usiamo per ora solo le informazioni relative al fuso orario locale.
PM10dataSA$Date.Local <- as.Date(PM10dataSA$Date.Local)Ora che i dati sono in formatO Date, R capisce che sono numeri salvati con un particolare formato e non stringhe.
In questo modo possiamo, ad esempio, vedere il periodo temporale che coprono usano i comandi min e max.
paste(min(PM10dataSA$Date.Local), which.min(PM10dataSA$Date.Local))[1] "2016-01-01 1"paste(max(PM10dataSA$Date.Local), which.max(PM10dataSA$Date.Local))[1] "2016-04-30 2559"Possiamo notare che il dataset non supera Aprile e che il max non è assunto nell’ultima osservazione, che sarebbe la 2867. Usando il comando tail() vediamo di capire perché.
tail(PM10dataSA$Date.Local, 25) [1] "2016-04-29" "2016-04-29" "2016-04-30" "2016-04-30" "2016-04-30"
[6] "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30"
[11] "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30"
[16] "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30"
[21] "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30" "2016-04-30"Fortunatamente vediamo che tutte le ultime 24 rilevazioni sono state fatte durante lo stesso giorno, pertanto il risultato di which.min() non è un errore. Ce lo potevamo aspettare, visto che le misurazioni sono ogni ora, ma è sempre meglio controllare.
Usando la funzione summary() estrapoliamo altre informazioni dal dataset.
summary(PM10dataSA) X.2 X.1 X State.Code County.Code
Min. : 1 Min. : 1 Min. :308088 Min. :6 Min. :9
1st Qu.: 646 1st Qu.: 771 1st Qu.:308858 1st Qu.:6 1st Qu.:9
Median :1291 Median :1442 Median :309529 Median :6 Median :9
Mean :1291 Mean :1464 Mean :309551 Mean :6 Mean :9
3rd Qu.:1936 3rd Qu.:2192 3rd Qu.:310279 3rd Qu.:6 3rd Qu.:9
Max. :2581 Max. :2867 Max. :310954 Max. :6 Max. :9
Site.Num Parameter.Code POC Latitude Longitude
Min. :1 Min. :81102 Min. :3 Min. :38.2 Min. :-120.7
1st Qu.:1 1st Qu.:81102 1st Qu.:3 1st Qu.:38.2 1st Qu.:-120.7
Median :1 Median :81102 Median :3 Median :38.2 Median :-120.7
Mean :1 Mean :81102 Mean :3 Mean :38.2 Mean :-120.7
3rd Qu.:1 3rd Qu.:81102 3rd Qu.:3 3rd Qu.:38.2 3rd Qu.:-120.7
Max. :1 Max. :81102 Max. :3 Max. :38.2 Max. :-120.7
Datum Parameter.Name Date.Local
WGS84:2581 PM10 Total 0-10um STP:2581 Min. :2016-01-01
1st Qu.:2016-02-02
Median :2016-03-01
Mean :2016-03-02
3rd Qu.:2016-04-02
Max. :2016-04-30
Time.Local Date.GMT Time.GMT Sample.Measurement
19:00 : 116 2016-01-02: 24 03:00 : 116 Min. : 1.00
20:00 : 116 2016-01-03: 24 04:00 : 116 1st Qu.: 4.00
06:00 : 113 2016-01-04: 24 02:00 : 113 Median : 7.00
07:00 : 113 2016-01-26: 24 05:00 : 113 Mean : 8.24
18:00 : 113 2016-01-28: 24 14:00 : 113 3rd Qu.:11.00
21:00 : 113 2016-01-29: 24 15:00 : 113 Max. :60.00
(Other):1897 (Other) :2437 (Other):1897
Units.of.Measure MDL Uncertainty
Micrograms/cubic meter (25 C):2581 Min. :4 Mode:logical
1st Qu.:4 NA's:2581
Median :4
Mean :4
3rd Qu.:4
Max. :4
Qualifier Method.Type Method.Code
Mode:logical FEM:2581 Min. :122
NA's:2581 1st Qu.:122
Median :122
Mean :122
3rd Qu.:122
Max. :122
Method.Name State.Name
INSTRUMENT MET ONE 4 MODELS - BETA ATTENUATION:2581 California:2581
County.Name Date.of.Last.Change DateTime.Local
Calaveras:2581 2016-05-06:619 2016-01-01 00:00:00: 1
2016-05-12:656 2016-01-01 01:00:00: 1
2016-06-22:618 2016-01-01 02:00:00: 1
2016-11-08:688 2016-01-01 03:00:00: 1
2016-01-01 04:00:00: 1
(Other) :2575
NA's : 1 Possiamo notare che alcune delle variabili sono trattate come numeriche anche se dovrebbero essere di tipo Factor. Possiamo o convertirle, o semplicemente tenerlo a mente qualora dovessimo lavorarci. Inoltre potremmo eliminare le colonne che non ci interessano usano l’assegnazione <- NULL. Ad esempio, dopo aver verificato che latitudine e longitudine sono le stesse in tutto il dataset, e corrispondono alla città di San Andreas, potremmo eliminarle. Va prima verificato che il dataset non contenga dati estranei, altrimenti eliminando una variabile, potremmo non essere più in grado di capirlo!
Tale verifica si può fare (in questo caso, ma può dipende dai dati) o usando il comando unique(), o anche leggendo i risultati di summary().
Volendo visualizzare i dati, che sono raccolti giornalmente, potrebbe fare comodo una colonna che riporta data e ora, in un formato che R riconosce. Possiamo ottenere tutto ciò con un solo comando. Va segnalato che il datset che è stato fornito, già contiene tale colonna, che è stata creata con il comando:
PM10dataSA$DateTime.Local <- as.POSIXct(paste(PM10dataSA[,c("Date.Local")] , PM10dataSA[,c("Time.Local")] ), format = "%Y-%m-%d %H:%M", tz="America/Los_Angeles",usetz=TRUE)Ora possiamo visualizzare i dati usando la funzione plot() aggiungendo delle opportune label lungo gli assi per migliorare la leggibilità.
plot(PM10dataSA$DateTime.Local, PM10dataSA$Sample.Measurement, xlab = "2016", ylab= "PM 10 (ug/m3)", cex = .5)
Poiché il grafico contiene un punto per ogni misurazione oraria, il grafico contiene moltissime informazioni.
Aggiungiamo una linea che indica il limite massimo giornaliero di PM 10 consentito (in Italia).
plot(PM10dataSA$DateTime.Local, PM10dataSA$Sample.Measurement, xlab = "2016", ylab= "PM 10 (ug/m3)", cex = .5)
abline(h = 50, col = "red")
# Per salviare l'immagine
#png(filename="../plot/DailyPM10.png")Si può notare che (fortunatamente) poche misurazioni superano i limiti consentiti.
2 Dati giornalieri
Vediamo di ridurle il numero di dati, aggregandoli per avere misurazioni giornaliere. Possiamo farlo usando la funzione aggregate() che ci permette di applicare una funzione ad un dataset indicando. Possiamo, ad esempio, decidere di salvare in un nuovo dataset le misure medie e massime per i dati di SA.
DailyVal <- aggregate(PM10dataSA$Sample.Measurement ~ Date.Local, data = PM10dataSA, FUN = mean)
DailyVal[,3] <- aggregate(PM10dataSA$Sample.Measurement ~ Date.Local, data = PM10dataSA, FUN = max)[2]
head(DailyVal)Possiamo migliorare la leggibilità del dataset cambiando i nomi delle variabili.
names(DailyVal)[2:3] <- c("mean daily PM10", "max daily PM 10")
head(DailyVal)Visualizziamo ora i dati medi e massimi giornalieri:
plot(DailyVal$Date.Local, DailyVal$`mean daily PM10`, xlab = "2016", ylab= "PM 10 (ug/m3)", cex = .8, main = "Concentrazione media giornaliera di PM10")
plot(DailyVal$Date.Local, DailyVal$`max daily PM 10`, xlab = "2016", ylab= "PM 10 (ug/m3)", cex = .8, main = "Concentrazione massima giornaliera di PM10")
abline(h = 50, col = "red")
Potrebbe essere interessante vedere se il giorno della settimana influenza la concentrazione di PM10, ad esempio in giorni con più traffico potremmo aspettarci più inquinamento.
typeof(DailyVal$Date.Local)[1] "double"#Sys.setlocale("LC_TIME","C") #Nomi
DailyVal$DayOfTheWeek <- weekdays(DailyVal$Date.Local, abbreviate = TRUE)
head(DailyVal)#ordiniamo i giorni in modo che la settimana inizi di lunedì
# ci sarà utile dopo
DailyVal$DayOfTheWeek <- ordered(DailyVal$DayOfTheWeek, levels=c( "Lun" , "Mar", "Mer" , "Gio", "Ven" , "Sab", "Dom" ))Vediamo di colorare il grafico precedente usando un colore diverso per ogni giorno della settimana.
plot(DailyVal$Date.Local, DailyVal$`max daily PM 10`, xlab = "2016", ylab= "PM 10 (ug/m3)", cex = .8, main = "Concentrazione massima giornaliera di PM10", col = DailyVal$DayOfTheWeek)
abline(h = 50, col = "red")
Vediamo che non sembra che i picchi siano raggiunti negli stessi giorni. Comunque una legenda andrebbe aggiunta. Vediamo con dei boxplot come la quantità media di polveri PM10 è distribuita rispetto ai giorni della settimana.
boxplot(`mean daily PM10` ~ DayOfTheWeek, data = DailyVal, ylab= "PM 10 (ug/m3)")
Aggiungere ad esempio una linea che indica le media delle misurazioni, aiuta a fare un confronto tra i dati.
boxplot( `mean daily PM10` ~ DayOfTheWeek, data = DailyVal, ylab= "PM 10 (ug/m3)")
abline(h = mean(DailyVal$`mean daily PM10`), col = "red")
Vediamo la data in cui è stato assunto il valore massimo:
DailyVal[DailyVal$`mean daily PM10` > 20,]Che sarà successo quel giorno? Potrebbe essere un errore di lettura?
Se volessimo visualizzare la distribuzione dei dati rispetto ai giorni, possiamo usare la funzione hist(). Purtroppo tale funzione non supporta la notazione ~, perciò dovremo specificare un grafico per giorno della settimana.
Dal grafico possiamo capire come si distribuiscono le misurazioni rispetto ai giorni della settimana. Si usi l’help per comprendere i parametri che sono stati usati.
2.1 Regressione (Esercizio)
Possiamo analizzare la relazione tra tempo e inquinamento usando un modello lineare? Proviamo e cerchiamo di capire se il modello è adeguato.
3 (Extra) Inquinamento in California
Decidiamo ora di concentrare la nostra attenzione sull’inquinamento in California e dopo aver letto il dataset che contiene i valori medi giornalieri di PM10, salviamo un sottoinsieme di dati che contiene solo dati relativi alla California.
PM10data <- read.csv("../data/daily_81102_2016.csv")
dim(PM10data)[1] 118885 29Vediamo che il file è molto grande, contiene molte righe e molte variabili. R impiega diverso tempo a leggerlo.
(PM10dataCAL <- PM10data[PM10data$State.Name == "California",])In questo caso il nome della variabile che ci interessa è Arithmetic.Mean.
Ora vogliamo visualizzare su una mappa quali sono lo città che superano la soglia di inquinamento consentita durante le osservazioni. In R, questa operazione è molto semplice usando i pacchetti maps o RgoogleMaps.
#install.packages("maps")
library(maps)
map("county", "california", xlim=c(-125,-114), ylim=c(32,43))
points( PM10dataCAL[PM10dataCAL$Arithmetic.Mean > 50,c("Longitude", "Latitude")] ,cex = .8, col = "red")
Se vogliamo usare una mappa di Goole Maps possiamo farlo usando il pacchetto dedicato.
library(RgoogleMaps)
map <- RgoogleMaps::GetMap(center="California", zoom = 6)
map1 <- plotmap( lat = PM10dataCAL[PM10dataCAL$Arithmetic.Mean > 50,c( "Latitude")], lon = PM10dataCAL[PM10dataCAL$Arithmetic.Mean > 50,c("Longitude")] ,map = map)
4 (Extra) Inquinamento a Los Angeles
Usiamo lo stesso dataset contenente i dati giornalieri relativi alle PM10 negli Stati Uniti e concentriamoci sulla sola città di Los Angeles.
PM10dataLA <- PM10data[PM10data$City.Name== "Los Angeles",]
dim(PM10dataLA)[1] 90 29Vediamo che sono state raccolte 90 osservazioni relative alla sola città di LA, cerchiamo di capire quando e con che frequenza sono state raccolte.
typeof(PM10dataLA)[1] "list"PM10dataLA$Date.Local[1:10] [1] 2016-01-01 2016-01-07 2016-01-13 2016-01-19 2016-01-25 2016-01-31
[7] 2016-02-06 2016-02-12 2016-02-18 2016-02-24
336 Levels: 2016-01-01 2016-01-02 2016-01-03 2016-01-04 ... 2016-12-02Vediamo che le date sono state salvate come liste. In realtà, come già visto, R prevede il formato Date che potrebbe essere utile e più maneggevole di una lista per fare “operazioni aritmetiche”. Convertiamo quindi la colonna in Date.
PM10dataLA$Date.Local <- as.Date(PM10dataLA$Date.Local)
typeof(PM10dataLA$Date.Local)[1] "double"La variabile ora risulta di tipo Date. Questo ci permette, ad esempio, di vedere quanti giorni passano tra una rilevazione e l’altra.
diff(PM10dataLA$Date.Local)Time differences in days
[1] 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[15] 6 6 6 6 6 6 7 5 6 6 6 6 6 6
[29] 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[43] 6 12 -270 6 6 6 6 6 6 6 6 6 6 6
[57] 6 6 12 6 6 6 6 6 6 6 6 6 6 6
[71] 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[85] 6 6 6 6 6Salta all’occhio che la frequenza delle registrazioni è circa ogni 6 giorni, ma compare un dato inatteso. Investighiamo meglio cosa è successo.
unique(PM10dataLA$Date.Local) [1] "2016-01-01" "2016-01-07" "2016-01-13" "2016-01-19" "2016-01-25"
[6] "2016-01-31" "2016-02-06" "2016-02-12" "2016-02-18" "2016-02-24"
[11] "2016-03-01" "2016-03-07" "2016-03-13" "2016-03-19" "2016-03-25"
[16] "2016-03-31" "2016-04-06" "2016-04-12" "2016-04-18" "2016-04-24"
[21] "2016-04-30" "2016-05-07" "2016-05-12" "2016-05-18" "2016-05-24"
[26] "2016-05-30" "2016-06-05" "2016-06-11" "2016-06-17" "2016-06-23"
[31] "2016-06-29" "2016-07-05" "2016-07-11" "2016-07-17" "2016-07-23"
[36] "2016-07-29" "2016-08-04" "2016-08-10" "2016-08-16" "2016-08-22"
[41] "2016-08-28" "2016-09-03" "2016-09-09" "2016-09-15" "2016-09-27"
[46] "2016-05-06" "2016-09-21"length(unique(PM10dataLA$Date.Local))[1] 47Capiamo che sebbene le osservazioni siano 90, non tutte si riferiscono a giorni differenti, quindi bisogna capire come gestire le misurazioni ripetute durante lo stesso giorno. Vediamo, ad esempio, cosa caratterizza le misurazioni ripetute del primo giorno dell’anno.
PM10dataLA[PM10dataLA$Date.Local== "2016-01-01",]Possiamo osservare che ci sono almeno due siti diversi dove vengono raccolti i dati. Questa osservazione segue sia la comparsa di due diversi Site.Num per la stessa data, sia le differenze nei valori di Latitude e Longitude. Vediamo se esistono solo due siti o più.
paste(unique(PM10dataLA$Site.Num) , unique(PM10dataLA$Latitude), unique(PM10dataLA$Longitude))[1] "1103 34.06659 -118.22688" "5005 33.9508 -118.43043" Abbiamo appurato che esistono due (e solo due) siti individuati da diverse coordinate geografiche.
Estrapoliamo ora alcune informazioni per le variabili numeriche usando il comando summary().
summary(PM10dataLA) State.Code County.Code Site.Num Parameter.Code POC
Min. :6 Min. :37 Min. :1103 Min. :81102 Min. :1.0
1st Qu.:6 1st Qu.:37 1st Qu.:1103 1st Qu.:81102 1st Qu.:1.0
Median :6 Median :37 Median :3054 Median :81102 Median :1.5
Mean :6 Mean :37 Mean :3054 Mean :81102 Mean :1.5
3rd Qu.:6 3rd Qu.:37 3rd Qu.:5005 3rd Qu.:81102 3rd Qu.:2.0
Max. :6 Max. :37 Max. :5005 Max. :81102 Max. :2.0
Latitude Longitude Datum Parameter.Name
Min. :33.95 Min. :-118.4 NAD83: 0 PM10 Total 0-10um STP:90
1st Qu.:33.95 1st Qu.:-118.4 WGS84:90
Median :34.01 Median :-118.3
Mean :34.01 Mean :-118.3
3rd Qu.:34.07 3rd Qu.:-118.2
Max. :34.07 Max. :-118.2
Sample.Duration Pollutant.Standard Date.Local
24 HOUR :90 PM10 24-hour 2006:90 Min. :2016-01-01
24-HR BLK AVG: 0 1st Qu.:2016-03-07
Median :2016-05-15
Mean :2016-05-14
3rd Qu.:2016-07-21
Max. :2016-09-27
Units.of.Measure Event.Type Observation.Count
Micrograms/cubic meter (25 C):90 Excluded: 0 Min. :1
Included: 0 1st Qu.:1
None :90 Median :1
Mean :1
3rd Qu.:1
Max. :1
Observation.Percent Arithmetic.Mean X1st.Max.Value X1st.Max.Hour
Min. :100 Min. : 6.00 Min. : 6.00 Min. :0
1st Qu.:100 1st Qu.:16.00 1st Qu.:16.00 1st Qu.:0
Median :100 Median :20.50 Median :20.50 Median :0
Mean :100 Mean :22.18 Mean :22.18 Mean :0
3rd Qu.:100 3rd Qu.:28.00 3rd Qu.:28.00 3rd Qu.:0
Max. :100 Max. :57.00 Max. :57.00 Max. :0
AQI Method.Code
Min. : 6.00 Min. :63
1st Qu.:15.00 1st Qu.:63
Median :19.00 Median :63
Mean :20.54 Mean :63
3rd Qu.:26.00 3rd Qu.:63
Max. :52.00 Max. :63
Method.Name
HI-VOL SA/GMW-1200 - GRAVIMETRIC :90
- : 0
BGI Inc. Model PQ200 PM10 - Gravimetric : 0
Hi Vol SSI Ecotech Model 3000 - Gravimetric: 0
HI-VOL-SA/GMW-321-B - GRAVIMETRIC : 0
HI-VOL-WEDDING-INLET - GRAVIMETRIC : 0
(Other) : 0
Local.Site.Name
LAX Hastings :45
Los Angeles-North Main Street :45
: 0
1 Site (Met Station and Hi-Vols BA1 and BA1-B) : 0
10 METER TOWER, AT NW CORNER, U OF I RESEARCH CENTER: 0
19th & Burt : 0
(Other) : 0
Address State.Name
1630 N MAIN ST, LOS ANGELES :45 California:90
7201 W. WESTCHESTER PARKWAY :45 Alabama : 0
151 NO SUNRISE BLVD, ROSEVILLE, CA : 0 Alaska : 0
201 ASHVILLE ROAD : 0 Arizona : 0
2220 NORTH STREET, ANDERSON, CA 96007: 0 Arkansas : 0
3337 Sandy Way, South Lake Tahoe, CA : 0 Colorado : 0
(Other) : 0 (Other) : 0
County.Name City.Name
Los Angeles:90 Los Angeles :90
Ada : 0 Aberdeen : 0
Adair : 0 Ajo : 0
Adams : 0 Ak-Chin Village: 0
Alamosa : 0 Alamosa : 0
Albany : 0 Albuquerque : 0
(Other) : 0 (Other) : 0
CBSA.Name Date.of.Last.Change
Los Angeles-Long Beach-Anaheim, CA:90 2016-12-20:90
: 0 2016-04-11: 0
Aberdeen, SD : 0 2016-04-15: 0
Albuquerque, NM : 0 2016-05-11: 0
Allentown-Bethlehem-Easton, PA-NJ : 0 2016-05-13: 0
Altoona, PA : 0 2016-05-18: 0
(Other) : 0 (Other) : 0 Possiamo notare che, anche in questo caso, alcune delle variabili sono trattate come numeriche anche se dovrebbero essere di tipo Factor. Possiamo o convertirle, o semplicemente tenerlo a mente qualora dovessimo lavorarci. Possiamo altrimenti rimuoverle come visto prima.
Per prima cosa, definiamo un campo di interesse. Se vogliamo, ad esempio, vedere l’andamento delle misurazioni medie di PM10 nei due siti, molte delle variabili non ci interessano. Possiamo investigare più nel dettaglio le variabili legate alla nostra indagine e non considerare le altre:
boxplot(Arithmetic.Mean ~ Site.Num, data = PM10dataLA)
Possiamo notare che uno dei due siti presenta delle misurazioni che in generale sono più alte.
Vediamo spazialmente dove sono posizionati i due punti di monitoraggio, questo potrebbe spiegare le differenze nelle misurazioni.
Dalla mappa sembra che uno dei due siti (cerchi neri) sia in centro mentre il secondo è più vicino alla spiaggia, dove possiamo aspettarci meno traffico.
Per avere una comprensione migliore visualizziamo una mappa differente, che mostra sia le strade che l’immagine satellitare della zona.
Vediamo che il secondo punto di osservazione, sebbene registri valori più bassi, probabilmente perché più vicino alla spiaggia, è comunque vicino all’aeroporto. Questo sicuramente influenza la qualità dell’aria e potrebbe spiegare perché sebbene più bassi, i valori non siano bassissimi. Non ci stupirebbe se altri parametri misurati fossero addirittura peggiori rispetto a quelli registrati in centro.
Senza opportuni strumenti di visualizzazione dei dati, non saremmo potuti giungere a queste conclusioni.